﻿using System;
using System.IO;
using System.Diagnostics;

namespace wicDIG
{
    /// <summary>
    /// Matrix 的摘要说明。
    /// 实现矩阵的基本运算
    /// </summary>
    public class Matrix
    {

        //构造方阵
        public Matrix(int row)
        {
            m_data = new double[row, row];

        }
        public Matrix(int row, int col)
        {
            m_data = new double[row, col];
        }
        public Matrix(double[,] d)
        {
            m_data = d;
        }
        //复制构造函数
        public Matrix(Matrix m)
        {
            int row = m.Row;
            int col = m.Col;
            m_data = new double[row, col];

            for (int i = 0; i < row; i++)
                for (int j = 0; j < col; j++)
                    m_data[i, j] = m[i, j];
        }

        /*
       //分配方阵的大小
        //对于已含有内存的矩阵，将清空数据
        public void SetSize(int row)
        {
            m_data = new double[row,row];
        }

        
        //分配矩阵的大小
        //对于已含有内存的矩阵，将清空数据
        public void SetSize(int row,int col)
        {
            m_data = new double[row,col];
        }
        */

        //unit matrix:设为单位阵
        public void SetUnit()
        {
            for (int i = 0; i < m_data.GetLength(0); i++)
                for (int j = 0; j < m_data.GetLength(1); j++)
                    m_data[i, j] = ((i == j) ? 1 : 0);
        }

        //设置元素值
        public void SetValue(double d)
        {
            for (int i = 0; i < m_data.GetLength(0); i++)
                for (int j = 0; j < m_data.GetLength(1); j++)
                    m_data[i, j] = d;
        }

        // Value extraction：返中行数
        public int Row
        {
            get
            {

                return m_data.GetLength(0);
            }
        }

        //返回列数
        public int Col
        {
            get
            {
                return m_data.GetLength(1);
            }
        }

        //重载索引
        //存取数据成员
        public double this[int row, int col]
        {
            get
            {
                return m_data[row, col];
            }
            set
            {
                m_data[row, col] = value;
            }
        }

        //primary change
        //　初等变换　对调两行：ri<-->rj
        public Matrix Exchange(int i, int j)
        {
            double temp;

            for (int k = 0; k < Col; k++)
            {
                temp = m_data[i, k];
                m_data[i, k] = m_data[j, k];
                m_data[j, k] = temp;
            }
            return this;
        }


        //初等变换　第index 行乘以mul
        Matrix Multiple(int index, double mul)
        {
            for (int j = 0; j < Col; j++)
            {
                m_data[index, j] *= mul;
            }
            return this;
        }


        //初等变换 第src行乘以mul加到第index行
        Matrix MultipleAdd(int index, int src, double mul)
        {
            for (int j = 0; j < Col; j++)
            {
                m_data[index, j] += m_data[src, j] * mul;
            }

            return this;
        }

        //transpose 转置
        public Matrix Transpose()
        {
            Matrix ret = new Matrix(Col, Row);

            for (int i = 0; i < Row; i++)
                for (int j = 0; j < Col; j++)
                {
                    ret[j, i] = m_data[i, j];
                }
            return ret;
        }

        //binary addition 矩阵加
        public static Matrix operator +(Matrix lhs, Matrix rhs)
        {
            if (lhs.Row != rhs.Row)    //异常
            {
                System.Exception e = new Exception("相加的两个矩阵的行数不等");
                throw e;
            }
            if (lhs.Col != rhs.Col)     //异常
            {
                System.Exception e = new Exception("相加的两个矩阵的列数不等");
                throw e;
            }

            int row = lhs.Row;
            int col = lhs.Col;
            Matrix ret = new Matrix(row, col);

            for (int i = 0; i < row; i++)
                for (int j = 0; j < col; j++)
                {
                    double d = lhs[i, j] + rhs[i, j];
                    ret[i, j] = d;
                }
            return ret;

        }

        //binary subtraction 矩阵减
        public static Matrix operator -(Matrix lhs, Matrix rhs)
        {
            if (lhs.Row != rhs.Row)    //异常
            {
                System.Exception e = new Exception("相减的两个矩阵的行数不等");
                throw e;
            }
            if (lhs.Col != rhs.Col)     //异常
            {
                System.Exception e = new Exception("相减的两个矩阵的列数不等");
                throw e;
            }

            int row = lhs.Row;
            int col = lhs.Col;
            Matrix ret = new Matrix(row, col);

            for (int i = 0; i < row; i++)
                for (int j = 0; j < col; j++)
                {
                    double d = lhs[i, j] - rhs[i, j];
                    ret[i, j] = d;
                }
            return ret;
        }


        //binary multiple 矩阵乘
        public static Matrix operator *(Matrix lhs, Matrix rhs)
        {
            if (lhs.Col != rhs.Row)    //异常
            {
                System.Exception e = new Exception("相乘的两个矩阵的行列数不匹配");
                throw e;
            }
            Matrix ret = new Matrix(lhs.Row, rhs.Col);
            double temp;
            for (int i = 0; i < lhs.Row; i++)
            {
                for (int j = 0; j < rhs.Col; j++)
                {
                    temp = 0;
                    for (int k = 0; k < lhs.Col; k++)
                    {
                        temp += lhs[i, k] * rhs[k, j];
                    }
                    ret[i, j] = temp;
                }
            }

            return ret;
        }


        //binary division 矩阵除
        public static Matrix operator /(Matrix lhs, Matrix rhs)
        {
            return lhs * rhs.Inverse();
        }

        //unary addition单目加
        public static Matrix operator +(Matrix m)
        {
            Matrix ret = new Matrix(m);
            return ret;
        }

        //unary subtraction 单目减
        public static Matrix operator -(Matrix m)
        {
            Matrix ret = new Matrix(m);
            for (int i = 0; i < ret.Row; i++)
                for (int j = 0; j < ret.Col; j++)
                {
                    ret[i, j] = -ret[i, j];
                }

            return ret;
        }

        //number multiple 数乘
        public static Matrix operator *(double d, Matrix m)
        {
            Matrix ret = new Matrix(m);
            for (int i = 0; i < ret.Row; i++)
                for (int j = 0; j < ret.Col; j++)
                    ret[i, j] *= d;

            return ret;
        }

        //number division 数除
        public static Matrix operator /(double d, Matrix m)
        {
            return d * m.Inverse();
        }

        //功能：返回列主元素的行号
        //参数：row为开始查找的行号
        //说明：在行号[row,Col)范围内查找第row列中绝对值最大的元素，返回所在行号
        int Pivot(int row)
        {
            int index = row;

            for (int i = row + 1; i < Row; i++)
            {
                if (m_data[i, row] > m_data[index, row])
                    index = i;
            }

            return index;
        }
        public static int Mod(int a, int m)
        {
            return (a - a / m * m) < 0 ? a - a / m * m + m : a - a / m * m;
        }
        /// <summary>
        /// return the submatrix with row1,col1 as the upperleft element and row2,col2 as the lowerright element
        /// if (row1,col1) == (row2,col2) then it will return a recombination of the 
        /// </summary>
        /// <param name="row1"></param>
        /// <param name="col1"></param>
        /// <param name="row2"></param>
        /// <param name="col2"></param>
        /// <returns></returns>
        public Matrix SubMat(int row1, int col1, int row2, int col2)
        {
            Matrix ret = new Matrix(Mod(row2-row1,Row)+1,Mod(col2-col1,Col)+1);
            for (int r = 0; r < ret.Row; r++)
                for (int c = 0; c < ret.Col; c++ )
                    ret[r, c] = this[Mod(r+row1,Row), Mod(c+col1,Col)];
            return ret;
        }
        /// <summary>
        /// Calculate Determinant of Square matrices
        /// </summary>
        /// <returns></returns>
        public double Det()
        {
            
            if (!IsSquare())
                return 0;
            else if (Row == 1)
                return this[0, 0];
            else
            {
                double result = 0;
                int sign = 1;
                for (int c = 0; c < Col; c++)
                {
                    if(this[0,c] != 0) //save time
                        result += this[0, c] * sign * SubMat(1, Mod(c + 1, Col), Row - 1, Mod(c - 1, Col)).Det();
                    sign *= -1;
                }
                return result;
            }
        }

        //inversion 逆阵：使用矩阵的初等变换，列主元素消去法
        /// <summary>
        /// this inverse algorithm is yet problematic
        /// </summary>
        /// <returns></returns>
        public Matrix Inverse()
        {
            if (Row != Col)    //异常,非方阵
            {
                System.Exception e = new Exception("求逆的矩阵不是方阵");
                throw e;
            }
            return new Matrix(ReverseMatrix(this.m_data, this.Row));
        }//end Inverse

        #region
        private double[,] ReverseMatrix( double[,] dMatrix, int Level )
        {
            double dMatrixValue = MatrixValue( dMatrix, Level );
            if( dMatrixValue == 0 ) return null;
                   
            double[,] dReverseMatrix = new double[Level,2*Level];
            double x, c;
            // Init Reverse matrix
            for( int i = 0; i < Level; i++ )
            {
                for( int j = 0; j < 2 * Level; j++ )
                {
                    if( j < Level )
                        dReverseMatrix[i,j] = dMatrix[i,j];
                    else
                        dReverseMatrix[i,j] = 0;
                }
         
                dReverseMatrix[i,Level + i ] = 1;
            }
         
            for( int i = 0, j = 0; i < Level && j < Level; i++, j++ )
            {
                if( dReverseMatrix[i,j] == 0 )
                {
                    int m = i;
                    for( ; dMatrix[m,j] == 0; m++ );
                    if( m == Level )
                        return null;
                    else
                    {
                        // Add i-row with m-row
                        for( int n = j; n < 2 * Level; n++ )
                            dReverseMatrix[i,n] += dReverseMatrix[m,n];
                    }
                }
         
                // Format the i-row with "1" start
                x = dReverseMatrix[i,j];
                if( x != 1 )
                {
                    for( int n = j; n < 2 * Level; n++ )
                        if( dReverseMatrix[i,n] != 0 )
                            dReverseMatrix[i,n] /= x;
                }
         
                // Set 0 to the current column in the rows after current row
                for( int s = Level - 1; s > i;s-- )
                {
                    x = dReverseMatrix[s,j];
                    for( int t = j; t < 2 * Level; t++ )
                        dReverseMatrix[s,t] -= ( dReverseMatrix[i,t]* x );
                }
            }
         
            // Format the first matrix into unit-matrix
            for( int i = Level - 2; i >= 0; i-- )
            {
                for( int j = i + 1 ; j < Level; j++ )
                    if( dReverseMatrix[i,j] != 0 )
                    {
                        c = dReverseMatrix[i,j];
                        for( int n = j; n < 2*Level; n++ )
                            dReverseMatrix[i,n] -= ( c * dReverseMatrix[j,n] );
                    }
            }
         
            double[,] dReturn = new double[Level, Level];
            for( int i = 0; i < Level; i++ )
                for( int j = 0; j < Level; j++ )
                    dReturn[i,j] = dReverseMatrix[i,j+Level];
            return dReturn;
        }
         
        private double MatrixValue( double[,] MatrixList, int Level )
        {
            double[,] dMatrix = new double[Level, Level];
            for( int i = 0; i < Level; i++ )
                for( int j = 0; j < Level; j++ )
                    dMatrix[i,j] = MatrixList[i,j];
            double c, x;
            int k = 1;
            for( int i = 0, j = 0; i < Level && j < Level; i++, j++ )
            {
                if( dMatrix[i,j] == 0 )
                {
                    int m = i;
                    for( ; dMatrix[m,j] == 0; m++ );
                    if( m == Level )
                        return 0;
                    else
                    {
                        // Row change between i-row and m-row
                        for( int n = j; n < Level; n++ )
                        {
                            c = dMatrix[i,n];
                            dMatrix[i,n] = dMatrix[m,n];
                            dMatrix[m,n] = c;
                        }
         
                        // Change value pre-value
                        k *= (-1);
                    }
                }
         
                // Set 0 to the current column in the rows after current row
                for( int s = Level - 1; s > i;s-- )
                {
                    x = dMatrix[s,j];
                    for( int t = j; t < Level; t++ )
                        dMatrix[s,t] -= dMatrix[i,t]* ( x/dMatrix[i,j] );
                }
            }
         
            double sn = 1;
            for( int i = 0; i < Level; i++ )
            {
                if( dMatrix[i,i] != 0 )
                    sn *= dMatrix[i,i];
                else
                    return 0;
            }
            return k*sn;
        }

        #endregion

        //determine if the matrix is square:方阵
        public bool IsSquare()
        {
            return Row == Col;
        }

        //determine if the matrix is symmetric对称阵
        public bool IsSymmetric()
        {

            if (Row != Col)
                return false;

            for (int i = 0; i < Row; i++)
                for (int j = i + 1; j < Col; j++)
                    if (m_data[i, j] != m_data[j, i])
                        return false;

            return true;
        }

        //一阶矩阵->实数
        public double ToDouble()
        {
            Trace.Assert(Row == 1 && Col == 1);

            return m_data[0, 0];
        }

        //convert to string
        public override string ToString()
        {

            string s = "";
            for (int i = 0; i < Row; i++)
            {
                for (int j = 0; j < Col; j++)
                    s += string.Format("{0} ", m_data[i, j]);

                s += "\r\n";
            }
            return s;

        }


        //私有数据成员
        private double[,] m_data;

    }//end class Matrix
}
